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Summary 

An analysis was made of the error in the min- 
imum altitude of a geometric ray from an orbiting 
spacecraft to the Sun. The dominant errors result 
from two sources: determining the position of the 
spacecraft and determining the position of the Sun. 

The sunrise and sunset errors are highly corre- 
lated and are opposite in sign. With the ephemeris 
generated for the SAGE I instrument data reduction, 
these errors can be as large as 200-350 m (lcr) after 
7 days of orbit propagation. The bulk of this er- 
ror results from errors in the position of the orbiting 
spacecraft rather than errors in computing the posi- 
tion of the Sim. These errors, in turn, result from the 
discontinuities in the ephemeris tapes resulting from 
the orbital determination process. Data taken from 
the end of the “definitive” ephemeris tape are used 
to generate the “predict” data for the time interval 
covered by the next arc of the orbit determination 
process. The predicted data are then updated by us- 
ing the tracking data. The growth of these errors is 
very nearly linear, with a slight nonlinearity caused 
by the beta angle. An approximate analytic method 
is given that predicts the magnitude of the errors and 
their growth in time with reasonable fidelity. 

Introduction 

A number of scientific reports have recently been 
cited in the press that describe sundry, sometimes 
conflicting, scenarios depicting the status and evolu- 
tion of the ozone content of the Earth’s atmosphere. 
Since ozone effectively filters out most of the po- 
tentially harmful ultraviolet radiation from the Sun, 
the sanctity and preservation of this most important 
trace gas, and its future in our atmosphere, is ob- 
viously of prime importance to all who inhabit this 
planet. 

In October of 1986, the National Aeronautics 
and Space Administration (NASA), in collaboration 
with the National Oceanic and Atmospheric Admin- 
istration (NOAA), the Federal Aviation Adminis- 
tration (FA A), the World Meteorological Organiza- 
tion (WMO), and the United Nations Environmen- 
tal Panel (UNEP) formed the Ozone Trends Panel. 
This panel, which consists of more than 100 scien- 
tists from all over the world, has been charged with 
the responsibility for evaluating the sundry findings, 
which involves the evaluation of ozone measurements 
from ground-based as well as satellite-borne instru- 
ments, and for determining which, if any, of the re- 
ported trends, catastrophic or otherwise, are in fact 
supportable under close scientific scrutiny and exam- 
ination from a common scientific base. 


As a result of the requirements of this panel, 
detailed error analyses were made of every facet 
involved in the generation of the final data prod- 
uct from each experiment (each of the analyses 
was originally reported independently in the tech- 
nical literature, and some of them led to the press 
reports mentioned earlier) — e.g., instrument mea- 
surement noise and calibration errors and their 
changes with time, data conversion, and data in- 
version algorithms — for both the long-established in- 
struments (e.g., Dobson spectrophotometers, rocket 
sounds) and the more recently developed satellite in- 
struments (SBUV, TOMS, LIMS, MAPS, SAGE I, 
and SAGE II). 

The SAGE I and SAGE II instruments are both 
used to perform solar occultation experiments in 
which the instrument scans across the solar disk 
at sunrise and sunset (with respect to the orbit- 
ing satellite) and measures the intensity of the so- 
lar radiation as it passes through various depths of 
the Earth’s atmosphere (see, e.g., McCormick 1982, 
and Mauldin et al. 1985). Scans outside the at- 
mosphere, exoatmospheric scans, provide the refer- 
ence radiation with which subsequent measurements 
are ratioed, and hence these instruments are essen- 
tially self-calibrating. Scans whose lines of sight pass 
through the Earth’s atmosphere have some of their 
solar energy removed, either by scattering (Rayleigh 
and aerosol) or by the absorption of some of this 
energy by molecules that are optically active in the 
bandpass of the instrument. The ratios of the atmo- 
spheric scans on a given pass to the exoatmospheric 
scans for the same pass provide a sequence of trans- 
mission functions that can be inverted to yield the 
altitude distribution of the density of the active ma- 
terial (vertical profiles), i.e., of ozone, in the present 
discussion. These vertical profiles can in turn be in- 
tegrated to yield the total vertical burden of ozone, 
for direct comparison with the Dobson data. 

Two of the several uncertainties in the data re- 
duction process that affect the accuracy of the end 
product generated by SAGE I and II are related to 
the ray that passes from the scan point on the Sun 
to the instrument — the line of sight, or tangent point 
ray — namely, (1) the minimum altitude of this ray 
above the surface of the oblate Earth and (2) the 
geographic location of this point. These are more 
fundamentally related to three other questions that 
arise frequently in orbital mechanics and/or spherical 
geometry problems: 

1. Where is the spacecraft located with respect 
to a pseudoinertial coordinate system whose origin is 
at the center of the Earth and whose X-axis coincides 
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with the vernal equinox (to be defined more precisely 
later)? 

2. Where is the center of the solar disk in the 
same coordinates? 

3. Where is the Greenwich meridian located in 
the same coordinate system — i.e., what is the Green- 
wich Sidereal Time? 

A first order analysis of the second set of ques- 
tions, and how these relate to the determination 
of the tangent point altitude, is the subject of the 
present paper. 

Error Analysis of Ephemeris Data 

The SAGE I and II ephemeris data are generated 
by the Goddard Space Flight Center (GSFC) by us- 
ing highly sophisticated trajectory computation and 
orbit determination programs. The spacecraft posi- 
tion and velocity vectors are generated in an inertial 
coordinate system defined as the mean of 1950.0 co- 
ordinates. These data are generated at 1-min in- 
tervals (plus other significant event times such as 
apogee and perigee passage, equatorial crossings). 
The ephemeris data are generated for a period of 
8 days, copied onto magnetic tape and sent to the 
users at the Langley Research Center (LaRC). 

This definitive ephemeris is the result of a number 
of processes, outlined in flowchart form in figure 1, 
each of which can inject errors or uncertainties into 
the final data product. 

In the first process, an orbital state vector at 
some time t — 0 is numerically integrated to produce 
8 days of ephemeris data at nominally l-min inter- 
vals. The mathematical model used in this step by 
the GSFC is quite sophisticated and includes many 
gravitational and nongravitational force terms whose 
magnitudes have been accurately determined by em- 
pirical methods from years of tracking Earth satel- 
lites. The resulting orbit is called a “predictive” or- 
bit and represents as good a set of ephemeris data 
as can be produced a priori by numerically solving 
the six-dimensional set of differential equations that 
contain the mathematical model of the force system 
acting on the satellite. 

From the predict orbits, a set of simulated track- 
ing data sets for each station on the tracking network 
is generated. The position of each tracking station 
is input as accurately as possible — its latitude, lon- 
gitude, and distance above or below the reference 
ellipsoid — and the spacecraft range, range rate, ele- 
vation and azimuth angles and their rates, and the 
spacecraft rise and set times relative to the station 
are accurately computed. 

The stations then proceed to track the actual 
spacecraft (see fig. 2). The real station tracking 


data are compared with the predicted data. The 
differences — called residuals — are then fed into a sta- 
tistical filter, usually a weighted-least-squares or a 
maximum likelihood filter, and corrections to the ini- 
tial conditions computed in such a way that the re- 
sulting computed (or modeled) orbit now minimizes 
some predetermined function of the residuals. This 
orbit is now called a “definitive” orbit, and it is the 
one whose ephemeris appears on the GSFC tapes. 
This represents the best set of orbital data available 
for the time period spanned by the ephemeris tapes. 
Because of the way the data are reduced and the 
ephemeris is computed, the orbit that this process 
yields is not continuous — there are discontinuities in 
the orbit state vector from one ephemeris tape to the 
next. This is shown schematically in figure 3. In the 
figure, the orbit is depicted as a sequence of disjointed 
arcs, while in reality each “arc” is about 105 orbits 
(7 days at approximately 15 orbits per day). 

There are three fundmental error sources that 
produce ephemeris data that deviate somewhat from 
reality. 

1. Imperfect mathematical models: The mathe- 
matical model used to generate the predict orbit is 
necessarily an approximation. Nature is a perfect 
analogue model, and any attempt to reproduce her 
results with a finite dimensional numerical model is 
necessarily going to culminate in only an approxima- 
tion to the real world situation. Hence, even if the 
initial conditions of a satellite orbit were perfectly 
known, the results of the model simulation would 
immediately begin to deviate from the real state of 
the satellite and the results would only be of accept- 
able accuracy until such time as the deviation became 
perceptible to some measurement process. This time 
might vary from a few tens of minutes (for a spherical 
Earth model) to several days or weeks for the GSFC 
models, depending on the mathematical model used 
and the purpose for which the ephemeris data are 
generated. 

2. Tracking errors: The tracking hardware used 
by the on- net tracking stations is subject to both me- 
chanical and electronic errors. These can only track 
to certain accuracies in range, range rate, and the an- 
gular data. The accuracies quoted in the literature 
for these parameters are a few meters (5-15) in range, 
6-30 mm/sec in range rate, 20 arc-sec in angle, and 
0.01 arc-sec/sec for angular rates (for angular rates 
less than about 500 arc-sec/sec). 

3. Orbit determination program errors: A com- 
plete orbit determination program would be very 
large, very complex, and very slow. There are many 
numerical approximations made in the production 
versions of these programs to speed up the solution 
process, and these necessarily reduce the accuracy 
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of the computations to varying degrees. Typical of 
the types of approximation made are numerical fi- 
nite difference methods to compute the many partial 
derivatives used in the inversion, and the lineariza- 
tion of some nonlinear effects. Other small errors are 
introduced in simplifying some of the mathematical 
models describing the physical processes occurring in 
the substructure components (e.g., simple methods of 
refraction and the effects of the Sun and the Moon 
on the satellite orbit). Finite word computers also 
introduce round-off errors which amplify in time. 

Determination of Errors in Initial Conditions 

The end result of all these errors and uncertainties 
manifests itself in the following way. The GSFC 
prepares the definitive orbit tapes and sends these 
to the working groups at LaRC. Each tape contains 
8 days of ephemeris data. The eighth day of one 
tape is an overlap of the first day of the next tape 
in the time sequence. The errors mentioned above 
show up when one compares the individual elements 
of the state vectors at the same times from these 
overlapping data sets. 

At the time the present study was initiated, there 
were 139 such overlaps on tapes containing SAGE I 
ephemeris data and 49 overlaps on the SAGE II 
tapes. The state vectors corresponding to the same 
data and time of day for each of the 139 SAGE I 
tapes were differenced, and the differences formed 
into a variance-covariance matrix according to the 
following standard maximum-likelihood estimator: 

cov(Ar) = E[{ r 8 - ri) (r 8 - r!) r ] 

= ^ “ r b)( r 8, " r b) T (!) 

t'=i 

where rg is the state vector at the beginning of 
day 8 of one tape, ri is the state vector at the 
corresponding time from day 1 of the next tape, 
and Ar = rg — ri. For SAGE I, N = 139, and 
for SAGE II, N = 49. This procedure resulted in 
the following set of standard deviations in the state 
vector : 



SAGE I 

SAGE II 

m 

326.0 

195.7 

III 

391.0 

172.5 

cr 2 , m 

281.4 

127.4 

ct-c, m/sec .... 

0.35 

0.18 

cry, m/sec .... 

0.32 

0.10 

(j Z) m/sec .... 

0.39 

0.21 


The reason for the apparent improvement in the 
SAGE II data is that in generating the 8-day tape for 


SAGE I a full 8 days of tracking data were processed 
and a definitive orbit was determined for the full 
time span. In the SAGE II processing, two separate 
4-day definitive orbits were computed independently 
and placed on one 8-day tape. From here on, only 
SAGE I results will be discussed, with the reasonable 
assumption that the SAGE II results will be about 
half those for SAGE I (see appendix A) . 

It is not possible to determine what proportion 
of these errors comes from errors in the state vector 
of day 7 of the first tape and how much comes from 
day 1 of the next tape. It is, therefore, assumed here 
that the actual errors at the beginning of day 1 of 
any tape are very small and that essentially all the 
error accumulated at day 7 of a given tape is a result 
of small errors at day 1 being propagated forward 
through 7 days of time. The only justification for 
this assumption, as will be shown, is that the error 
growth at the end of a 7-day period is rather close to 
the magnitude of the errors observed in the above- 
mentioned comparison of the state vector. Note that 
7 days of propagation time is used, as t = 0 of a 
given tape is equivalent to the beginning of day 8 of 
the previous tape, which occurs at the end of the 7th 
day plus one orbit. 

The next step is to determine what values of 
the standard deviations in the uncertainties of the 
initial conditions are required at t — 0 such that 
the above values of uncertainties are realized after 
7 days of error propagation. One of the SAGE I tapes 
was picked arbitrarily, and the initial conditions at 
t = 0 of day 1 were assumed to represent nominal 
conditions. An initial set of standard deviations was 
assumed and a random number generator used to 
determine 20 sets of initial conditions for a set of 
perturbed orbits. Each of these were run out for the 
full 7 days, and each hour the differences between the 
perturbed and nominal orbits were collected. The 
means and standard deviations were computed at 
each hour and the results plotted. The initial guess 
at the standard deviations in the initial conditions 
was adjusted until the standard deviations of the 
state vector elements after 7 days agreed reasonably 
well with the table above, computed from the real 
ephemeris data. These standard deviations at t = 0 
were found to be 

a Xo — 0.3227 m 

a yo - 0.3870 m 

a Zo - 0.2781 m 

a x 0 = 0.3465 x 10 -3 m/sec 

(j y o = 0.3168 x 10~ 3 m/sec 

o Zo = 0.3861 x 10 -3 m/sec 
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Figure 4 shows a typical plot of the resulting 
history of the x-component of the state vector errors. 
It is seen that the midpoint of the oscillation at 
t — 168 hr (7 days) is very close to the 326 m called 
for (see also appendix B). 

The next step, using the standard deviations 
at t = 0 just found, was to generate a number 
of perturbed orbits and compute statistics on the 
pertinent parameters. A sample size of 100 was 
chosen here and the parameters whose statistics were 
generated were 

1. The difference between the tangent point alti- 
tudes at sunrise for the perturbed and nominal or- 
bits. The nominal altitude chosen here was 10 km. 
A few check runs showed that the altitude difference 
was essentially independent of the magnitude of the 
nominal altitude. 

2. The same as above for the sunset events. 

3. The difference between the sunrise tangent 
point altitude and the very next sunset tangent point 
altitude on the same orbit. 

The rise and set altitudes were determined as 
follows. The initial conditions were taken from the 
SAGE I ephemeris tapes. These initial conditions 
were assumed to be a nominal set for the beginning 
of each 7-day period. For a nominal tangent point 
altitude of 10 km, the time of every rise and set event 
during the ensuing 7-day period was computed from 
a program based on Buglia (1986) and written to a 
data file, along with the (3 angle (see fig. 5) at each 
event time. The nominal initial conditions were then 
fed into a second program which, using a random 
number generator, computed 100 sets of random 
initial conditions, Gaussian distributed with zero 
mean and standard deviation given above for t = 0. 
A third program then computed the tangent point 
altitudes for each rise and set event for each orbit for 
the same times as determined for the nominal orbit 
(see fig. 6) . For each of the 100 perturbed orbits, the 
three statistics described above were determined for 
each event and the results plotted. 

The entire process was repeated for 25 consecu- 
tive weeks in order to get a reasonable variation of the 
(3 angle into the results. As shown in figure 5, (3 is the 
angle between the direction to the Sun and the or- 
bital plane, measured positive in the direction of the 
angular momentum vector of the orbit. For (3 = 90° 
(a situation not encountered on either SAGE I or 
SAGE II, the maximum useful value of (3 for the 
SAGE instruments is about 66°), the orbital plane 
is normal to the Sun direction, and it is apparent 
that even a large error along the orbit will have little 
effect on the altitude of the tangent point altitude. 
For /? = 0°, on the other hand, the Sun “sees” the 


orbit edge-on, and here small errors along the orbit 
path will result in large errors in tangent point al- 
titude. Figure 5 shows this /3 = 0° situation more 
clearly, and this is the geometry that would produce 
the maximum error in tangent point altitude. 

A few histograms were plotted showing the dis- 
tribution of the quantity (ft/j — hjyji), where hji is 
the rise altitude of the perturbed orbit and hj^g — 
10 km, the nominal rise altitude. Two of these 
are presented in figures 7 and 8, showing that with 
only a little imagination these appear to be normally 
distributed with essentially zero mean. Two cases 
with extreme variation were chosen here for display 
purposes. The distributions are thereby assumed 
normal, and hence it makes sense from here on to 
characterize their distributions by use of standard 
deviations. 

The standard deviation histories of the three 
quantities mentioned earlier are shown from a typical 
run on the next three figures, figures 9-11. These are 
the standard deviations in the quantities (hR-h^n), 
( hg — Jins), and (hft — kg), where hg is the set 
altitude and hpjg — 10 km, the nominal set altitude. 
The statistic (Hr — hg) does not include any nominal 
data (except implicitly in the time, of course) but is 
the difference between the rise altitude and the very 
next set altitude in the same perturbed orbit. 

These curves are typical of the histories of the 
noted statistical quantities for small to moderate 
variations in the (3 angle (the effects of f3 will be 
shown later). In these figures, f3 increases from -33° 
at t = 0 to about -4° at t = 168 hr. The standard 
deviation increases from some very small — but not 
zero — value at t = 0 to something of the order of 
340 m for both (h,R — h^ft) and (hg — hpfg) (figs. 9 
and 10) and to about twice this (as expected) for the 
quantity (hg — hg) (fig. 11), and this increase is very 
nearly linear (see appendix A). 

Much of the effect of errors in the initial coor- 
dinates manifests itself as a very small error in the 
mean angular velocity of the spacecraft and conse- 
quently shows up as a displacement error along the 
flight path in the nominal orbit plane. As shown in 
appendix A, for zero (3 angle, the angle between the 
nominal and perturbed radius vectors at time t is 
given approximately by 


2 \a 3 / a 


( 2 ) 


where a is the semimajor axis of the nominal orbit 
and fi is the gravitational constant for the Earth, 
398600.64 km 3 /sec 2 . For example, for the orbits 
of figures 9-11, a = 6963 km and for a moderate 
Aa = 0.001 km, the above equation gives an angle 
9 after 7 days of 1.415 x 10 — 4 rad. This yields in 
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a down- track error of just under 1 km. From the 
relation (appendix A) 

Ah = - (r SCyN cos 4> n ) 6 (3) 

one computes a rise and set altitude error of about 
390 m, which is somewhat higher than that shown in 
figures 9-10 but is the correct order of magnitude and 
is very close in magnitude to some of the standard de- 
viations computed for other time periods. Figure 12 
shows that for nearly circular orbits, a “fast” tim- 
ing error would cause an increase in the sunrise alti- 
tude and a corresponding decrease, of about the same 
magnitude, in the sunset altitude. This is borne out 
in the figures and is also obviously the reason why 
the error in the difference (/i# — hg) is about twice 
the error in either sunrise or sunset altitude. This 
also indicates that the rise and set altitude differ- 
ence should be highly correlated from one orbit to 
the next, and this is amply illustrated in figure 13, 
which shows a time history of (h# — hg) for one of 
the 100 perturbed orbits used in the statistical data 
compilation. For this particular choice, this quantity 
is negative over the entire period and the high degree 
of correlation between orbits is shown clearly — given 
a few hours history of this difference, the difference 
for the next few days is quite predictable. It is im- 
portant to emphasize this fact that, starting with a 
given ephemeris tape from the GSFC, a bias will de- 
velop between the calculation of the rise and set alti- 
tudes and that this bias will increase approximately 
linearly over the 7-8 days of computation time. For 
the example shown, this bias grew to about 700 m in 
7 days, the computed set altitude being higher than 
the computed rise altitude (at the nominal times), 
indicating a fast time error (the semimajor axis of 
the computed orbit is slightly smaller than the (un- 
known) semimajor axis of the real orbit). 

The variation of the altitude error with the j3 
angle is shown in figure 14. The standard deviation 
of the quantity cr(hji — hg) is shown in the upper half 
of the figure and that for a{hji — h^n) is shown in 
the lower half, both plotted against /?, at t = 168 hr 
for each of the 25 weeks of data. To get the point 
at 0 = 0°, the mean of the quantity cr(Ah)/cos /? 
was computed using the results at t — 168 hr. The 
solid curves are then plots at A cos f3, where A is the 
amplitude for each of the statistics described above. 

Solar Ephemeris Errors 

A number of procedures were used in the SAGE I 
and II data reduction algorithms to accurately 
compute the position of the center of the solar disk 
with respect to the center of the Earth and subse- 
quently with respect to the orbiting spacecraft. The 


coordinate system chosen for all subsequent calcula- 
tions is the mean of date system, in which the X- 
axis (the vernal equinox) is determined by the inter- 
section of the Equator with the ecliptic, precession 
being taken into account (nutation is neglected). The 
scan plane of the SAGE instruments contains the 
spacecraft and the centers of the Earth and the Sun. 
The positions of the top and the bottom of the Sun 
as seen in the scan plane are computed, with due 
allowance made for the refraction of the ray as it 
passes through the atmosphere. The algorithms ul- 
timately developed are completely self-sufficient in 
that they do not need any external ephemerides 
or other input — the only inputs needed are a cal- 
endar date and a Greenwich Mean Time, and the 
corresponding position and velocity vectors of the 
spacecraft. 

All the equations and methods used in these 
algorithms, except the aberration corrections, are 
outlined and discussed in Buglia (1988). 

Definition of Minimum Altitude 

Assume a quasi-inertial coordinate system with 
the origin at the center of the Earth (fig. 15). The 
Z-axis is along the rotational axis of the Earth, and 
the equatorial plane passes through the origin and is 
normal to the Z-axis. The X-axis, or vernal equinox, 
is located by definition as the point on the equatorial 
plane where the geocentric orbit of the Sun crosses 
the Equator from south to north. This occurs about 
21 March of each year. The Y-axis lies also in the 
equatorial plane in such a way as to complete a right- 
hand system. This is called a quasi-inertial system 
because the vernal equinox is not fixed in inertial 
space, but moves slowly with time (the precession of 
the equinox — see, e.g., Buglia 1988). 

The mathematical shape describing the surface 
of the Earth (sphere, oblate spheroid, geoid) is given 
in this coordinate system. Also given are three po- 
sition coordinates and three velocity coordinates of 
an orbiting spacecraft at some specific time, and the 
three position coordinates of some celestial body (in 
this study, the center of the disk of the Sun), both 
relative to the Earth-centered inertial system. From 
these parameters one can define a vector pointing 
from the spacecraft to the celestial body. If one starts 
at the spacecraft and makes one’s way along this vec- 
tor, or ray, and looks vertically downward (in the 
plane formed by the centers of the three objects — 
the centers of the Earth, Sun, and spacecraft), one 
finds that a minimum altitude above the Earth’s sur- 
face occurs somewhere along this ray. The point on 
the ray at which this minimum occurs is called the 
tangent point , and the point on the surface of the 
Earth immediately below the tangent point is called, 
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appropriately enough, the subtangent point. The dis- 
tance between the tangent point and the subtangent 
point is the altitude of the tangent point and is the 
minimum altitude referred to in the remainder of the 
present paper. The computation of the minimum 
height is trivial for a spherical Earth, but gets a bit 
more complicated when made for an oblate spheroid 
(see, e.g., Buglia 1988). A subroutine in the SAGE I 
and II data reduction software computes the magni- 
tude of this minimum distance and the inertial coor- 
dinates of the subtangent point. The geographic co- 
ordinates (ordinary geodetic latitude and longitude) 
are also computed. 


Error Sources in Minimum Height Calculation 

It can be seen from the above definition of the 
tangent point that there are three major geometric 
sources that contribute to errors in determining the 
minimum height of the observation ray above the 
Earth’s surface and two additional corrections that 
have to be made to yield the position of the Sun as 
seen at the spacecraft. The geometric error sources 
are 

1. The shape assumed for the surface of the 
Earth. 

2. Uncertainties in the geocentric position of the 
spacecraft. 

3. Uncertainties in the geocentric position of the 
center of the Sun. 

The additional corrections that need to be made are 

4. Planetary aberration due to the motion of the 
Earth around the Sun. 

5. Orbital aberration due to the motion of the 
spacecraft around the Earth. 

6. Atmospheric refraction. Uncertainties due to 
refraction uncertainties were not included in the 
present paper. See Chu 1983. 

1. Shape of the Earth. For a point on the surface 
of the rotating Earth, the potential of gravity, W, is 
the sum of the potential of the gravity force (due to 
the mass distribution of the Earth) and the centrifu- 

gal potential (due to the rotation of the Earth) (see, 

for example, Heiskanen and Moritz 1967, ch. 2) 

, oo n 
n=lm=0 

x [Jnm cos (mA) + K nm sin(mA)]| 

+ ^lr 1 2 zos 2 0 (4) 


in which 

p Earth gravity constant, 

398600.64 km 3 /sec 2 * 

r geocentric radius to a point on the 

surface, km 

a e equatorial radius of the Earth, 

6378.1600 km 

9 geocentric latitude, deg 

A longitude from Greenwich, deg 

u e angular velocity of the Earth’s rota- 

tion, 7.292 115855 x 10 -5 rad/sec 
where P nrn are the standard associated Legendre 
polynomials, and J nm and K nm are integrals involv- 
ing the mass distribution of the Earth. These are 
the gravitational constants mentioned earlier. These 
have been empirically determined by the GSFC from 
satellite tracking. 

For the SAGE I and II software, the assumption 
was made that the potential of gravity was rotation- 
ally symmetric about the north polar axis. This re- 
moves the longitudinal dependence from equation (4) 
and allows it to be written in the simpler form involv- 
ing only zonal harmonics 



r 


OO 

E(y)' JnPn(sm9) 



( 5 ) 


This dynamical form was used in some of the SAGE 
software, with the first six zonal harmonics being 

Ji = 0 

J 2 = 1082.6271 x HT 6 
J 3 = -2.5358868 x 10~ 6 
J 4 = -1.6246180 x 10 -6 
J 5 = -0.22698599 x 10 -6 
J 6 = 0.54518572 x 10~ 6 

The radius of the constant potential surface repre- 
sented by equation (5) was computed for the Earth’s 
surface potential (from eq. (5) at the Equator) and 
compared with the radius computed from a simple 
ellipsoid with equatorial radius of 6378.160 km and 
polar radius of 6356.775 km. These are shown in 
comparison in the table on the following page. It 
is seen that the maximum difference occurs at 45° 
latitude and is less than 82 m. This was considered 
small, and in all the algorithm segments involving 
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0, deg 

r e (ellipsoid), km 

r w (potential), km 

r w - r e , m 

0 

6378.14000 

6378.14000 

0 

5 

6377.97906 

6377.97485 

-4.21 

10 

6377.50100 

6377.48831 

-12.69 

15 

6376.72000 

6376.69553 

-24.47 

20 

6375.65925 

6375.62110 

-38.15 

25 

6374.35033 

6374.29822 

-52.11 

30 

6372.83230 

6372.76760 

-64.69 

35 

6371.15063 

6371.07617 

-74.46 

40 

6369.35588 

6369.27554 

-80.33 

45 

6367.50223 

6367.42048 

-81.75 

50 

6365.64588 

6365.56719 

-78.69 

55 

6363.84337 

6363.77170 

-71.67 

60 

6362.14980 

6362.08819 

-61.61 

65 

6360.61717 

6360.56740 

-49.77 

70 

6359.29271 

6359.25517 

-37.54 

75 

6358.21737 

6358.19105 

-26.31 

80 

6357.42447 

6357.40715 

-17.31 

85 

6356.93864 

6356.92714 

-11.50 

90 

6356.77500 

6356.76550 

-9.50 


the geometric Earth surface the ellipsoid of revolu- 
tion was used in place of the surface defined by the 
harmonic potential. The quoted values of the ellip- 
soidal radii, the J-coefficients, and the gravitational 
constant are the latest (as of March 1986) received 
from the GSFC. These data are for the Northern 
Hemisphere only. There is no north-south symme- 
try because of the nonvanishing of the odd-J terms 
(the J3 term, for example, gives the “pear shape” 
to the Earth). The difference in Earth radius at the 
corresponding north-south latitudes may be as high 
as 30-35 m, but in no case does the potential surface 
radius described from equation (5) deviate by more 
than 80 m or so from that of the ellipsoid. 

2. Errors in the position of the orbiting spacecraft. 

The spacecraft ephemeris computations are made in 
the data reduction software package by numerically 
integrating the equations of motion of the spacecraft 
in a quasi-inertial coordinate system fixed at the 
center of the Earth. The gravity field used in the 
model is the sixth-order zonal model described above. 

Ephemeris data, deduced from actual spacecraft 
tracking data, are received from the GSFC on mag- 
netic tapes as discussed earlier. These data are given 
at nominally 1-min intervals in the mean 1950.0 coor- 
dinate system and span an 8-day time interval. The 
adjective “mean” here is the standard astronomical 
meaning and indicates that the precessional effect of 
the Earth’s rotation has been included in the coor- 
dinate computation, but the nutational effects have 


not. The initial conditions for the ephemeris integra- 
tion are extracted from this tape for the time cor- 
responding to the nearest minute of time preceding 
the beginning of the data taking event. The initial 
state vector is transformed from the mean 1950.0 
coordinates to the mean of date coordinates, the 
date as specified on the GSFC tape, using the stan- 
dard transformation equations of spherical astron- 
omy (see, for example, Taff 1985). The equations of 
motion are then integrated in this coordinate system 
at 60 equal time intervals throughout the data taking 
period, which generally ranges from 70-90 sec. 

Near a graze event (0 near 66°) the tangent 
height decreases toward zero as the Sun begins to 
set, just touches the horizon, and immediately begins 
to increase as the Sun rises. These events may last 
as long as 10-12 min. Even for this length of time, 
the root sum square differences between the position 
and velocity vectors as computed from the zonal- 
only model described above and the GSFC results, 
accepting the GSFC results as a standard, are only 
of the order of 10 m and 20 mm/sec, respectively. 
These differences are deemed quite acceptable. 

3. Solar ephemeris calculations.. The computa- 
tion of the geocentric position of the center of the 
Sun for any day of any year, and any time of day, is 
accomplished by the use of highly accurate series ex- 
pressions that give the mean longitude, mean anom- 
aly, and true anomaly of the Sun, and the eccentricity 
of the Earth orbit, all computed with respect to the 
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mean equator and equinox of date (see algorithms in 
Buglia 1988). These allow the computation of the 
right ascension and declination of the center of the 
Sun with respect to the mean equinox and equator 
of date. 

In a similar way, the Mean Greenwich Sidereal 
Time is computed from formulas that give the MGST 
at 0 hr GMT for any date and are then corrected for 
time of day. 

4 . Planetary aberration. There is a slight shift in 
the position of the Sun as seen at the Earth’s center 
because of the finite speed of light. It takes about 
8.3 min for photons from the Sun to reach the Earth, 
so we really see the Sun where it was 8 min earlier. 
The ratio of the Earth’s orbital speed (30 km/sec) to 
the speed of light (300000 km/sec) is about 0.0001, 
and since the Earth’s velocity is essentially normal 
to the Earth-Sun line, the maximum displacement 
in the solar position due to planetary aberration is 
of the order of 20 arc-sec. Most of this shows up as 
an essentially constant decrease of 20 arc-sec or so in 
the right ascension of the Sun, with a periodic shift 
in declination of amplitude 9 arc-sec and period one 
year. This correction is made in the SAGE I and II 
software by using the classical aberration correction 
formulas of spherical astronomy (e.g., Smart 1977). 
The correction as applied here is the reverse of that 
usually made in classical astronomy. There, the cor- 
rection is applied to the “observed” coordinates to 
yield the “true” position. Here, the “true” coordi- 
nates are calculated, and the correction is added to 
yield the “observed” position. 

5. Orbital aberration. In addition to the Earth 
orbital motion about the Sun, the spacecraft orbiting 
about the Earth introduces an additional aberration 
effect due to its own velocity. To be strictly correct, 
the velocities of the spacecraft with respect to the 
Earth and that of the Earth with respect to the Sun 
should be added vectorially, and the component of 
this resultant normal to the line of sight from the 
spacecraft to the Sun used to compute the total 
aberration effect, which is then apportioned among 
the right ascension and declination in accordance 
with the classical formulas. 

In the SAGE I and II software, the orbital aber- 
ration has been neglected. The orbital velocity of 
the spacecraft is approximately 7.5 km/sec. This is 
about one-fourth of the orbital velocity of the Earth, 
and hence the maximum aberration shift caused by 
the orbital motion would be about 5 arc-sec. How- 
ever, this maximum value would only be realized if 
the spacecraft were traveling in a direction normal 
to the Earth-Sun line. In fact, however, during pe- 
riods of observation of sunrise and/or sunset, the 


spacecraft is traveling essentially in a direction to- 
ward or away from the Sun, a condition that min- 
imizes the aberration shift. Even during the graze 
events, when the plane of the spacecraft is as nearly 
normal to the Earth-Sun line as it gets during ob- 
servations, the displacement due to aberration is of 
the order of 2 arc-sec. The other extreme case occurs 
when the orbital plane contains the Earth-Sun line 
and the Sun is on the Earth horizon. The angle be- 
tween the velocity vector and the Earth-Sun line is 
about 24°, and this produces an aberration shift of 
about 3 arc-sec. Now, 1 arc-sec angular displacement 
at the spacecraft produces an altitude error of about 
14 m, and hence the maximum altitude error due to 
neglect of their orbital aberration is only about 40 m. 

Accuracy of Computations of Coordinates of 
Sun Center 

To check the accuracy of these computations, a 
number of calculations were made with the software 
and compared with appropriate data in the Astro- 
nomical Almanacs for a number of years. 

The Mean Greenwich Sidereal Times and the ap- 
parent right ascension and declination of the center 
of the Sun were extracted from the pertinent Astro- 
nomical Almanacs for 0 hr GMT for the first of each 
month for the years 1979-1987 inclusive. 

The Greenwich Sidereal Time calculations were 
found to be very close. The computed mean side- 
real times differed from the published values by 
only 0.001-0.002 sec for the years 1979-1983. From 
1984 to 1987 the difference between the computed 
and published times jumped to 0.070 sec and re- 
mained fixed for these years. The reason for this 
change is that in 1984 the International Astronom- 
ical Union adopted, among other resolutions, the 
following changes in their procedures for computing 
ephemerides: 

1. A new standard equinox (Jan 1.5, 2000 instead 
of Jan 0.5, 1900). 

2. A more accurate and self-consistent set of 
planetary masses. 

3. Small corrections to the precession/nutation 
constants. 

4. General relativity is now included in the equa- 
tions of motion of the planetary orbits. 

The new constants are now included in the com- 
putation of the SAGE ephemeris parameters for the 
years 1984+, and the differences between the com- 
puted and published Greenwich Sidereal Times are 
now again down in the 0.001-0.002 sec range. 

The differences between the computed and pub- 
lished coordinates of apparent right ascension and 
declination are shown graphically in figures 16 and 17, 
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respectively. These data are plotted for the first day 
of each month for the years 1979-1987, inclusive. 

The mean (average) error in the right ascension is 
5.4 arc-sec, with a standard deviation of 9.6 arc-sec. 
The corresponding quantities for the declination are 
-0.7 arc-sec and 5.0 arc-sec. Since both of these un- 
certainties are transcendentally related to uncertain- 
ties in the tangent point height, it is perhaps more 
useful — and certainly more extreme — to consider the 
root sum square of these errors and assume that this 
result acts directly in the direction of maximum tan- 
gent point altitude error. Figure 18 displays this 
parameter for the same conditions as stated ear- 
lier. The mean uncertainty here is 10.3 arc-sec with 
a standard deviation of 8.5 arc-sec. Even if all this 
error acted in the scan direction (which, as stated 
above, it does not), the resultant uncertainty in the 
tangent point altitude would be 147 ± 116 m (1 <t). 

Thus, in an absolutely worst case scenario, if 
all of the above uncertainties acted together in the 
worst way, the mean uncertainty in the tangent point 
height due to uncertainties in the computation of the 
solar position would be 147 + 116 + 40 or 303 m. 
Considering that the errors attain their maxima at 
different times and that they act in different direc- 
tions, it is reasonable to conclude that the most likely 
uncertainties due to ephemeris error are considerably 
less than 303 m, perhaps a factor of 5-10 less. 

Conclusions 

The results of the present study of ephemeris er- 
rors contained on the SAGE I ephemeris tapes ob- 
tained from the GSFC permit the following conclu- 
sions to be drawn: 


1. The accuracy of the ephemeris tapes is at 
worst several hundred meters, and it is doubtful that 
much improvement could be realized by the GSFC. 

2. The standard deviation in the error between 
the true (but unknown) rise and set tangent point al- 
titudes and those computed from the GSFC ephemeris 
tapes for SAGE I can amount to as much as 200- 
350 m at the end of 168 hr of error propagation. 

3. The rise and set tangent point errors are gener- 
ally of opposite sign and equal in magnitude (for the 
nearly circular orbits of SAGE I), so that the stan- 
dard deviations in the difference between the rise and 
set tangent point heights is about twice that for the 
individual errors, or about 400-700 m. 

4. On a given orbit, the rise and set tangent 
point altitude errors are always of about the same 
order of magnitude (for the nearly circular orbits of 
SAGE I) and are always opposite in sign. The sign 
persists, at least over the 7 days studied here, and 
the magnitudes vary essentially linearly with time. 
There is, therefore, a very high degree of correlation 
between these altitudes over the time span covered 
by a single ephemeris tape. 

5. The uncertainties introduced in the computa- 
tion of the tangent point altitudes due to uncertain- 
ties in determining the position of the spacecraft and 
the position of the Sun in the appropriate coordi- 
nate system are considerably smaller than those due 
to the spacecraft ephemeris computation errors. A 
worst case scenario, totally at variance with physi- 
cal reality, indicates that this error should be smaller 
than about 300 m, with the real world error proba- 
bly being one-half to one order of magnitude smaller 
than this. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
December 8, 1988 
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Figure 1. Schematic of the process leading to the generation of SAGE I and II ephemeris tapes 









— — Actual orbit (unknown) 

Predict 

Definitive (best guess) 

Figure 2. Illustration of the three orbits: actual, predict, and definitive. 



Figure 3. Discontinuities between weekly ephemerides. 









Figure 7. Histogram of the errors in sunrise tangent point altitude at 168 hr using the ephemeris tape for 
the week of 10/28/79. 


Number 
of errors 



Figure 8. Histogram of the errors in sunrise tangent point altitude at 168 hr using the ephemeris tape for 
the week of 9/28/80. 




Figure 9. History of the standard deviation in the sunrise altitude error for the week of 11/25/79, The 
magnitude of the /? angle decreases from 33° to 4° during the 7-day period. 



0 12 24 36 48 60 72 84 96 108 120 132 144 156 .168 

Time, hr 


Figure 10. Same as figure 9 for the sunset altitude error. 





Figure 11. Same as figure 9 for the difference between sunrise and sunset altitudes on the same orbit. 



Ah R > 0 
Nominal rise 


Nominal set 
Ah s < 0 


Figure 12. Sketch showing the increase in sunrise and decrease in sunset tangent point altitudes for a 
“fast” time error. 
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Figure 14. Effect of the f3 angle on the standard deviation in (1) the difference between the perturbed 
and nominal sunrise and sunset tangent point altitudes, and on (2) the difference between rise and set 
altitudes. 






Figure 16. Difference between computed and published right ascension of the Sun, first of each month for 
the 9-year period 1979-1987 inclusive, arc-sec. 



Figure 17. Difference between computed and published declination of the Sun, first of each month for the 
9-year period 1979-1987, inclusive, arc-sec. 



Figure 18. Plot of the quantity [(Ara)^ + (A dec)^}^^ for first day of each month for the 9-year period 
1979-1987, inclusive, arc-sec. 
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Appendix A 

Derivation of an Approximate Expression 
for the Standard Deviation in the Tangent 
Point Altitudes 

If we assume nearly circular orbits and assume 
that f3 — 0°, then it is possible to derive an approx- 
imate expression for the standard deviation in the 
quantity (Hr — h^n) or (hg — hjyg) f rom figure 6(a). 
We will here call this quantity Ah, since to the de- 
gree of approximation used here, it is the same for 
both rise and set tangent point altitudes. From the 
figure 


Ah = r sc , p sin <p p - r SCtN sin <p N (Al) 

where r sc is the distance of the spacecraft from the 
Earth’s center, and 4> is defined as the angle from the 
local vertical to the center of the Sun. The angle 9 
(fig. 6(a)) is the angular distance between the nom- 
inal radius vector and the perturbed at the time at 
which the nominal vector reaches the position to pro- 
duce the nominal value of the tangent point height. 
(Subscripts N and p refer to nominal and perturbed 
parameters, respectively.) If At is the time from the 
beginning of the computation period (zero time at 
the beginning of the tape), then approximately 

9 « 9 0 + (rip - n;v)At (A2) 

where 0 o is the angular separation between the nom- 
inal and perturbed radius vectors at t = 0, and n is 
the mean angular rate in the orbit, 


6 0 will generally be a very small angle (of the order 
of a few hundredths of an arc-sec) and hence will be 
neglected in what follows. 

From figure 5, 


4>n — 4>p + 9 


Since the initial condition errors are very small, and 
the errors in the magnitude of the radius vector are 
very small (see fig. 15), then using equations (A3) 
and (Al) and assuming that 6 « 1 we get 

Ah = - 9 (r SC:N cos <f>) (A4) 

The nominal value of n in equation (A2) is 




N , 


1/2 


where p = 398600.64 km 3 /sec 2 , the gravitational 
constant of the Earth, and a is the semimajor axis of 
the nominal orbit. If we expand n p in a Taylor series 
about riiy, 


n p 


n N + 



Aa -f- . . . 


From equation (A3) 

dn 3 n;v 

da 2 a/y 

and equation (A2) becomes 


9 = 



At 


(A5) 


(A6) 


(A7) 


In terms of the magnitudes of the initial radius and 
velocity vectors, 


J_ _ 2 
a N ~ r , 


o M 


(A8) 


and hence 


which gives 


da ( as \ 

dr 0 V r 0 ) 

da _ 2 a2 N v ° 

dv a p 


2 'k 


(A9) 


da da 

Aa = r—A r 0 + -r—Av 0 
or 0 ov 0 

2 2 
a\j 0'%jV 0 A 

— 2—ifAr 0 + 2^ — A Vo 


(A10) 


Finally, with 


r<> = (zii + vl + zl) 


1/2 


A r 0 = w^-Ax 0 + —-A y 0 + |^A z 0 

OXo ^Vo vZo 


= — (x 0 Ax 0 -I- y 0 Ay 0 + z 0 Az 0 ) (All) 
r 0 


and similarly 


A v 0 = — [i 0 Ai 0 + y 0 Ay 0 + z 0 Az 0 ] (A12) 

Vo 
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n N = 


In equations (A4) and (A7), the only random 
variables are 9 and Aa, respectively. Hence, if we 



assume Gaussian distributions with zero means 


when approximately 


cr(Ah) = |r SCjJV cos <^v|cr(0) 


and 


0(6) = 


t 

2 


ct(u) 


(A13) 

(A14) 


in which from equation (A 10), assuming no corre- 
lation between initial position and initial velocity 
errors, 


° ^ \o 2 (xo) + a 2 {y 0 ) + 0 2 {z o ) | 

o = | ^ [cr 2 (i 0 ) + a 2 {y 0 ) + 0 2 (z o ) 


1/2 
1/2 

For SAGE I, we have nominally a = 6980 km and 


a(a) = 


2a 


N \ _2 


O ' ( r °) 


+ 


2(l]yVo 


1/2 


0 2 (v o ) 


(A15) 


and 


CT ( r °) = ~ {[x 0 cr(x 0 )} 2 + [y 0 o{yo)? 
+ [zoCr(z 0 )] 2 j 1 

p{v 0 ) = {[i 0 a(i 0 )] 2 + [y 0 <r{y 0 )Y 

1/2 


+ [ir 0 cr(ir 0 )] 2 1 
If we make the assumption that 

o(x 0 ) = <r(t/ 0 ) = a(z 0 ) = a 


(A16) 


(A17) 


with a similar expression for the velocity components, 
then equations (A16) and (A17) become simply 

o{r 0 ) = cr 
a(v Q ) = a 


a — 3.32 x 10 4 km 
a = 3.51 x 10 -7 km/sec 


which gives 


er(a) = 9.23 X 10 4 km 
a(9) = 1.3 X 10~ 7 rad 
<j(A/i) = 369 m 


From equations (A13) and (A14) we see that, to the 
order of approximation introduced in this analysis, 
cr(A h) is a linear function of time, a conclusion that is 
approximately borne out in figures 9-11. The slight 
departure from nonlinearity in these figures is due to 
the variation in the (3 angle during the 7 days of error 
propagation. Identical plots for other time periods 
show that cr(Ah) versus time is slightly concave 
upward when (3 decreases over the time period and 
slightly concave downward when (3 increases, showing 
that as (3 gets large the magnitude of the errors 
decreases, and vice versa, as heuristically developed 
in the text (see fig. 6). 
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Appendix B 

Application of / and g Series to Error 
Propagation 

The / and g series method can readily be used 
to develop analytic expressions for the propagation 
of errors along an orbit. These relations state that 
(see, e.g., Escobal 1965) 

r(t, t 0 ) = /(£, t 0 )r 0 + g(t, t 0 )v 0 (Bl) 


where r(<) is the position vector at any time t, given 
in terms of the initial position vector r 0 and the 
initial velocity vector v 0 . The / and g series are, 
in general, infinite series derived by expanding the 
radius vector in a Taylor series, replacing r by the 
expression for the acceleration due to gravity as soon 
as r appears, and collecting the coefficients of r a 
(/ series) and r 0 (g series). For the simple two- 
body problem, where the gravitational force of the 
attracting body is a central force, the / and g series 
can be expressed in exact closed form (Escobal 1965) 

/(Mo) = 1 — - [1 — cos(£ - £„)] (B2) 

r 0 

g(t, t 0 ) = (t-t 0 ) ~-[{E - E 0 ) - sin(£ - E 0 )\ (B3) 
n 

in which a is the semimajor axis of the orbit, E 0 
and E are the initial and current values of eccentric 
anomaly, and n is the mean angular rate in the orbit 



in which \i is the gravitational constant of the central 
body (for the Earth, y, = 398600.64 km/sec). 

For circular orbits, the / and g expressions reduce 
considerably, and, for example, the x-component of 
equation (Bl) can be written as 

x(t, t 0 ) = x 0 cos n(t — t a ) H — - sin n(t — t 0 ) (B5) 

n 

with similar expressions for y and z. These can in 
turn be differentiated to yield expressions for the 
velocity components. Only the x-component will be 
developed here. 

Suppose now that uncertainties are introduced 
into the initial conditions — Ax 0 , A y 0 , . . . , Az 0 . 
Equation (A10), with equations (All) and (A12), 
still defines Aa. Then 


dx(t,t 0 ) . dx(t, t 0 ) . . 
A x(t,t 0 ) = — ^ — -Ax 0 + — ^ — -Ax 0 

oXq dxo 


+ 


dx(t,t 0 ) dn ^ 


dn da 


(B6) 


With the help of equation (B5) and a little algebra, 
we can write equation (B6) as 


A x(t, t 0 ) = Asin(n At) + B cos (n At) (B7) 


where 


A = 


/ Ax 0 3 x 0 A a\ 
\ n 2 n a ) 



(B8) 


3 A ci. 

B = Ax 0 — -x 0 — -At (B9) 

2 CL 

and A a follows from the equations of appendix A 
already cited. 

If we define 

tani/> - — (BIO) 

A 

then we can write equation (B7) in the more compact 
form 


A x(t, t a ) = V A 2 + B 2 sin(n At + ip) (Bll) 

The variance a^.{t,t 0 ) does not assume such a 
simple form, and it is easier to evaluate it from the 
defining equation 



dx(t,t 0 ) 

dx 0 



2 


+ 


(dx(t,t a ) \ 2 


+ 


( dx(t, t 0 ) dn 
V dn da 



2 


(B12) 


Equations (Bll) and (B12) were evaluated using 
the same initial conditions as were used to generate 
figure 13, along with the variances cited in the main 
text. The results from these two equations are plot- 
ted in figures Bl and B2, respectively. The similarity 
between figure 4 of the main text and figure B2 can- 
not be missed. The small differences are due to an in- 
consistency in the use of equations (Bll) and (B12). 
The nominal Cartesian coordinates used in the def- 
initions of the A and B parameters, equations (B8) 
and (B9), were the same as those used in generat- 
ing the single orbit whose partial results are shown 
in figure 13. These nominal Cartesian coordinates 
do not produce quite a circular orbit, as assumed in 
the development of equations (Bll) and (B12), but 
rather an orbit with an eccentricity of 0.00807. Con- 
sequently, r a ± a, and v 0 is not the circular orbital 
velocity at radius a. These inconsistencies introduce 
small timing errors in the computation of the mean 
angular rate, which are mainly responsible for the 
differences between figures 4 and B2. 
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Figure Bl. Analytic computation of the deviation in the ^-component of spacecraft position due to initial 
perturbations equal to the standard deviations given in the text. Computed from equation (Bll). 
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Figure B2. Same data as figure 4, but computed using the same inputs as in figure Bl, and using 
equation (B12). 
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